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1; 
-TITLE MTHSSQRTR2 #9 Point Square Root routine 
- IDENT /1-015/ ; : MTHSQRTR2.MAR EDIT: RNH1015 
MARBAARAAAALAALALALAASAAAAAAEAALEALESALELA ALES ESAS ELE R REESE SERRE SESE SEES EEE SS | 
COPYRIGHT (c) 1978, 1980, 1982, 1984 


° BY 
DIGITAL at ety _GORPORATION, MAYNARD, MASSACHUSETTS. 
ALL RIGHTS RESER 
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°® % 
°® ® 
*® x 
** * 
*® ® 
°® 0 ® 
;* N ACCORDANCE WITH THE TERMS OF SUCH ante a AND WITH THE * 
* ION OF THE ABOVE COPYRIGHT NOTICE. THIS SOFTWARE OR ANY OTHER * 
;* THEREOF a NOT BE PROVIDED OR OTHERWISE MADE AVAILABLE TO ANY * 
3% PERSON. NO TITLE TO AND OWNERSHIP OF THE SOFTWARE IS HEREBY * 
;* TRANSFERRED. * 
°* * 
;* THE INFORMATION IN THIS SOFTWARE IS SUBJECT TO CHANGE WITHOUT NOTICE * 
;* AND SHOULD NOT BE CONSTRUED AS A COMMITMENT BY DIGITAL EQUIPMENT * 
;* CORPORATION. * 
°® © 
+8 oa 
°* * 
of 4 
°® ® 
2 * 


DIGITAL ASSUMES NO RESPONS 
SOFTWARE ON EQUIPMENT WHIC 
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; g FACILITY: MATH LIBRARY 
: "ABSTRACT: 
: : MTHSSORT _R2 is a special routine which is the same as MTH$SQRT except 


a faster non-standard JSB call is used with the argument in RO and no 
; registers are saved. 


S HISTORY: ; 
; AUTHOR: 
: Peter Yuo, 15-Oct-76: Version 01 


; MODIFIED BY: 


; Ol-1 Peter Yuo, 22-May-77 
; 01-2 Peter Yuo, 31-May-77 


| 
| 
| 
; VERSION: 01 | 
| 
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-SBTTL HISTORY ; Detailed Current Edit History 


3; Flo 
HISTO 


ALGORITHMIC DIFFERENCES FROM FP=11/C ROUTINE: none 
Edit History for Version 01 of MTH$SQRT 


te Code saving after code review 

1-2 ROTL shift in garbag e into pi onest bit. Use ASHL instead. 
ADDL netryctiop matt er ADJUST has gocSt 16 shanged pote ADDW to prevent 
overflow if R 6> = PFFP ang, 

\-3 Finish error "Lendl ite 10-June-1 
MTHSSERROR eoenere te — STGNAL. 

MTH$_... change 
Changed error Seedl ine "sathitien. Put error result in RO before 
wk ng MTHSSSIGNAL in order to allow user modify error result. 


OOOO oQooQooQooooQooooooooooooo <r 
SOOCOQOOCOCOOCOCOOOCOOOOOOOOooOoOooOo 
SOOGCOCOCOCCOCOCOCOOCOOOOOOOOOOOoOO -+--35 


SOoOOoCoooooooooooooooooooooo we 


NAME WIN CO OODNAUE UP 9 OONAUE Wry 


NSN NENA AA AAAAAA AMNION 


gi-6 Return -0.0 on negative arg. TNH 0-Dec-77 
01-7 Edit in Rich Lary’s code bums. JSB routine is now _R2. JMT 19-Jan-78 
01-9 Move .ENTRY symbol to module header. TNH 14-Aug-78" 
1-010 = Put version number >. sreneurs format: three iets edit 
numbers. Also ones the Mah ht notice BS 16-NOV-78 
1-011 - cheng, io SauRO0 cs to *RTHSK UROONEG das 07-DEC-78 
eH ~ o the T dire tive. 
1-015 - Sects pine Bt SBL 17- may-1979 
1-014 = Move MTHS$SQRT_R2 to separate module. JAW ge-tegttr¥.. 
1-015 = Change external references to G*. RNH 06-Oct-8 
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-SBTTL DECLARATIONS ; Declarative Part of Module 


: INCLUDE FILES: 


EXTERNAL SYMBOLS: 


-DSABL GBL 
eEXTRN = MTHSK_ SQUROONEG 
-EXTRN MTHSSSIGNAL 


EQUATED SYMBOLS: 
MACROS: none 
PSECT DECLARATIONS: 


-PSECT _MTHSCODE PIC,SHR,LONG, EXE ,NOWRT 
; program section for math routines 


OWN STORAGE: none 


3 § ~OOODQOGDOOOCO0OO000000 0009 CR c7 CoCo CDCDODCO 
MEW O OONAU EWN 3 O ODNAUE WN O ODNAULS UM 00 


1 
1 

0000 1 

0000 1 

0000 =«(1 CONSTANTS: 

0000 «(1 

0000 106; 

44 \ 3; Constants A and B chosen for k = odd 
13CDS5FD4 0900 109 ° LF_ODD_A_E63_ = “X13CD5FD4 
3€4A2018 8208 : LF-ODD_B_EM63 = *x3C4A2018 

4 ! : Constants A and B chosen for k = even 
F61A4015 6000 114° LF_EVEN_A = “XF61A4015 
4B8231FD7 0000 1 LF_EVEN_B_EM64 = *X4B231F D7 
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~SBTTL MTHSSQRTR2 = Standard Single Precision Floating SQRT 


p++ 
: FUNCTIONAL DESCRIPTION: 


= single precision floating point function 


SQRT(X) is computed using the following approximation technique: 


If X <= 0, error. Let X = iXi. 
Let X = 2e*K * F where F is the fractional part. 


If K = even, X = 2*#(2P) * F 
SQRT(X) = 2eeP * SORT(F), 1/2 =< F< 1 


If K = odd, X = 28(2P4+1) ® F = 208(2P+2) © (F/2) 
SQRT(X) = 2e#(P41) ® SQRT(F/2), 1/4 =< F/B < 1/2. 


Let F' = AtF + B, 


A = 0.453730314(octal), 
B = 0.327226214(octal), for K = even. 

= At(F/2) + B, 
A = 0.650117146(octal), 

4 B = 0.230170444(octal), for K = odd 
an 
K' =P, for K = even 
=P? 1 for K = odd. 


Let YO = 2**kK' * F* as a staight line spgren inet ten wthin the 
given interval using coefficients A and B which minimize the 
absolute error at the midpoint and endpoint. 

Starting with YO, two Newton-Raphson iterations are performed. 
YCn+1) = (1/2) * ( YEn) + X/YEn]) 


The relative error is < 10**-8, 


; CALLING SEQUENCE: 


sart.wf.v = MTH$SQRTR2(x.rf.r) 


; INPUT PARAMETERS: 


LONG = 4 ; define Longuerg multiplier 
x = 1 * LONG ; Contents of x is the argument 


IMPLICIT INPUTS: none 
OUTPUT PARAMETERS: 


VALUE: floating square root of the argument 


IMPLICIT OUTPUTS: none 
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1-015 1984 11 


ATHSSQRTRE = Standard Single Precision 6=SEP- 

COMPLETION CODES: none 

SIDE EFFECTS: 

Signoiss MTHS_SQUROONEG if xX < 0. 4 with reserved BET anAL) in RO (copied to 
@ signal mechanism vector CHF$L_MCH_RO/R1 by LIBS$S 


Associated message is: ‘SQUARE ROOT OF NEGATE E VALUE". Result is reserved 
operand -0.0 unless a user supplied (or any) error handler changes CHF SL _MCH_ RO/RI | 


735 ett on Macro V04-00 
715 CMTHRTL.SR CIMT MTHSQRTR2.MAR; 1 


oOo 


NOTE: This procedure disables floating point underflow, enables integer 
overflow, causes no floating overflow or other arithmetic traps, and 
preserves enables across the call. 
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u 4-00 
2 = JSB SQRT routine 6-SEP-1984 11:27:1 MTHRTL.SRCJMTHSQRTR2.MAR; 1 
-SBTTL MTHSSQRT_R2 = JSB SORT routine 
JSB SQRT = used by the standard, and directly. 


CALLING SEQUENCE: 
gove anything in RO:R2 


owe : input in RO 
JSB MTH$SQRT_R2 
return with result in RO 


=z 
“— 
z=" 
An-—- 
“oOo 
ow 
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| 
| 
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| 


Note: This routine is written to avoid any integer overflows, floating overflows, | 

floating underflows or divide by 0 conditions, whether enabled or not. | 
REGISTERS USED: 

RO = Floating argument then resuit | 

Ri - X saved for use during iteration | 

| 


Sete Ge Se Se Ge Ge Ge Ge Ge Ge Ge Ge Ge Se 


R2 = scratch 
JSB routine for SQRT 


ee 
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MOVE 
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19 
198 
194 
Hf 
199 
198 
8 199 
00 
8 01 
03 
0000 $04 
0 05 
0 86 
8 0 
08 
0000 510 
> sem F 11 RO, R1 ; test sign of X and save it in R1. 
| 55 15 0 1g BLEQ  ZERO_NEG : branch to ZERO_NEG if X =< 0 
0 14; x%>0 
0 2 : 
0005 216 Pos: 
7E D& 0005 17 CLRL -(SP) ; make room for 2nd half of 
0007 18 ; double Length operand 
52 50 3¢ 0007 19 MOVZWL RO, R2 : isolate low 16 bits ain exp,>fract) in R 
2s 94 QO0A 220 CLRB R ; R2 new has sion and left ? exp bits 
50 65 AA 000C 221 BICW Ree RO ; clear sign and left 7 exp bits 
50 895 QO0F ss¢ TSTB R ; check low bit of exp 
10 18 0011 22 BGEQ EVEN 3; and branch if 1 
50 13CD5FD4 8F 44 0013 24 MULF #LF_ODD_A_E63, RO ; add 64 (half of bias) to (exponent-2) | 
OO1A $ : and start approximation calc 
50 3C€4A2018 8F 40 OO1A 5 ADDF § #LF_ODD_B_EM63, RO : RO = (first approx) * 2**-64 
13. «11 + 1 ssi BRB ADJOST 3 go adjust 
00 ; 229 EVEN: 4 
50 2000 8F Ad 00 230 ADDW #*x2000, RO : exp is 0 - make it 64 (2**-64) for legaiit | 
50 F61A4015 BF 44 0028 31 MULF #LF_EVEN_A, RO 
50 4B231FD7 BF 40 3p f 3 — ADDF #LF_EVEN_B_EM64, RO ; RO = (first approx) * 2**-64 | 
52 52 1F 9C 0036 3 ROTL #31, R2, R2 ; divide R2 (exp+bias) by 2, 
i 2 : giving (exp/2+64) 
50 52 ad A 6 ADDW R2, RO ; insert exp/2 in first approx and 
Bgep ¢ 3; re-bias it. 
00% 33 : first iteration - single precision is sufficient 
52 51 50 47 D 41 DIVF3 RO, R1, R2 ; R2 = 
50 52 40 0041 s ADDF 2k ; RO = YO + x/Y0 
50 0080 8F A2 0044 4 SUBW #*x80, RO ; RO = Y1_= (1/2)(¥O + X/¥0) 
984 re 3 no overflow possible 
b49 rt) 3; second iteration, do in double precision to get truncated( rather than 
049 47 ; rounded) result. 
049 48 ; 
} 
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51 ge 30 
55 3 
50 0080 8F 
FD 
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7 00° 8F 
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q Point Square Root =O ke ab eb 9 91: 
2 = JSB SQRT routine 6-SEP=-1984 

49 33: CLRL R 

29 $36 Divod rs, R1 

5 

38 

54 PUSHL Ri 

55 CLRL R1 

2$ DIVD3 RO, (SP)+, R1 

5 DDF R1, RO 

28 SUBW #*x80, RO 

59 SQRTX: RSB 

60 

61; xX =< 0 

$6 : 

63 ZERO_NEG: 

64 BEQL SQRTX 

65 PUSHL (SP) 

66 MOVZBL #MTHS$K_SQUROONEG, -(SP) 
67 ASHL #15, #T, RO 
268 

269 
$9 
si CALLS #2, G*MTHSSSIGNAL 

r RSB 

74 

75 ~END 
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ee 
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SR 


ro V04-00 
MTHRTL.SRCIMTHSQRTR2.MAR; 1 


lower part (xX) = 0 
g ide Y1 into X with low-order 
bits of Y1 garbage. This doesn't 
effect accuracy, since Y1 innacurate 
anyway. 
convert x and place on stack 
clear low part of Y1 
Ror"! nigher 

= 4 er_pa 
RO = SQr 2s (T/2 
return, ®) = result 


35 
15 
Ow 
iv 
2 

ff 


return with RO = result = 0 

return PC from JSB routine 
condition value 

RO = result = reserved operand -0.0 


; RO goes to signet mechanism vector 
(CHFSL_M 0 


CH_RO/R1) so error handler 
can modify the result. 


; ro error and use real user's PC 
in 


ependent of CALL vs JS 


B 
; return = RO restored from CHFS$L_MCH_RO/R1 
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| MTHSSQRTR2 ; Floating Point Square Root routine obs Se 9 01:51:35 
Symbol table 6-SEP=-1984 11:27:15 
ADJUST 89000 Bp 4 
EVEN 0000 aR 1 
F _LEVEN_A = F6IAG 

LF_EVEN B EM64 = sGesiF0 
LF-ODD_K_€63_ = 13CDS5FD4 
LF ODD_8_EM63 = TaD 
L = 00000004 
MTHS$$S1IGNAL rereeret X 00 
MTHSK_SQUROONEG ‘*tteneeee x 00 
MTHSSORT_R2 00800000 RG 01 
POS 00000005 R 01 
SQRTX eile TH R 01 
ZERO_NEG OOOO05SA R 01 

$eereneaenec wen arsoe + 

! Psect synopsis ! 

peeaer ener wroesees + 
PSECT name Allocation PSECT No. Attributes 
- ABS, 00000000 <( 0.) 00 ¢ 0.) NOPIC USR CON ABS LCL NOSHR NOEXE NO 
_MTHSCODE QOOOO06E ( 110.) O1¢ 1.) PIC USR CON REL LCL SHR_- EXE 

henennaneemnnnanadenmenns 
! Performance indicators ! 

Phase Page faults CPU Time Elapsed Time 
Initialization 3 00:00:00.07 00:00:01.65 
Command processing 12 00:00:00.69 00:00:03.66 
Pass 1 79 00:00:00.78 00:00:04.39 
Symbol table sort 0 00:00:00.00 00:00:00.00 
Pass 2 61 00:00:00.71 00:00:02.94 
Symbol table output ¢ 00:00:00.02 00:00:00.02 
Psect synopsis output 00:00:00.02 00:00:00.02 
Cross-reference output 0 00:00:00.00 00:00:00.00 
Assembler run totals 306 00:00:02.29 00:00:12.68 


The working set Limit was 900 pages. 

3664 bytes (8 pages) of virtual memory were used to buffer the intermediate code. 

There were 10 pages of symbol table space allocated to hold 14 non-local and 0 local symbols. 
5 source lines were read in Pass 1, producing 8 object records in Pass 2. 

1 page of virtual memory was used to define 1 macro. 


eee re em ee meme eee ee eee een ea} 


Macro Library name 


~$255$DUA28: CSYSLIBISTARLET.MLB;2 0 
0 GETS were required to define 0 macros. 


There were no errors, warnings or information messages. 
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| 
| 
| 
RD NOWRT NOVEC BYTE 
| 
| 
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) 


RD NOWRT NOVEC LONG 
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MACRO/ENABLE=SUPPRESSION/DI SABLE=(GLOBAL , TRACEBACK) /LIS=LISS:MTHSQRTR2/OBJ=OBJ$:MTHSQRTR2 MSRCS$:MTHJACKET/UPDATE=(ENHS:MTHJACKET) +MS 


ORPORATION 
OPRIETARY 


AH-BT13A-SE DI 
VAX/VMS V4.0. CO 


PMENT 
AND P 


80263 


244 AH BT13A-SE : ENT CORPORATION 
VAX/VMS V4.0 ND PROPRIETARY 


